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ABSTRACT 

(N ■ 

We discuss the dissipation of turbulent kinetic energy E^ in the global interstellar medium 
["•». ' (ISM) by means of two-dimensional, MHD, non- isothermal simulations in the presence of model 

radiative heating and cooling. We argue that dissipation in two dimensions is representative of 
that in three dimensions as long as it is dominated by shocks rather than by a turbulent cascade. 
Contrary to previous treatments of dissipation in the ISM, in this work we consider realistic, 
stellar-like forcing: energy is injected at a few isolated sites in space, over relatively small scales, 
and over short time periods. This leads to the coexistence of forced and decaying regimes in 
the same flow, to a net propagation of turbulent kinetic energy from the injection sites to the 
decaying regions, and to different characteristic dissipation rates and times in the forced sites 
and in the global flow. 

We find that the ISM-like flow dissipates its turbulent energy rapidly. In simulations with 
forcing, the input parameters are the radius of the forcing region, the total kinetic energy 
ek each source deposits into the flow, and the rate of formation of those regions, Sob- The 
global dissipation time id depends mainly on If. We find that for most of our simulations id 
is well described by a combination of parameters of the forcing and global parameters of the 
flow: id w itrms/ (^k/)j where u rms is the rms velocity dispersion, is the specific power of each 
forcing region, and / is the filling factor of all these regions. In terms of measurable properties 
of the ISM, id ^ (S g }«^ ms /(ekSoB), where (S g ) is the average gas surface density; for the solar 
neighborhood, id > 1-5 x 10 7 yr. The global dissipation time is consistently smaller than the 
crossing time of the largest energy-containing scales, suggesting that the local dissipation time 
near the sources must be significantly smaller than what would be estimated from large-scale 
quantities alone. 

In decaying simulations, we find that the kinetic energy decreases with time as E^(t) oc t~ a , 
where a « 0.8-0.9. This result can be translated into a decay with distance £ when applied to 
the mixed forced+decaying case, giving E^ cx g- 2a /( 2 - a ) a t large distances from the sources. 

Our results, if applicable in the direction perpendicular to galactic disks, support models of 
galaxy evolution in which stellar energy injection provides significant support for the gas disk 
thickness, but do not support models in which this energy injection is supposed to reheat an 
intra-halo medium at distances of up to 10-20 times the optical galaxy size, as the dissipation 
occurs on distances comparable to the disk height. However, this conclusion is not definitive until 
the effects of stratification on our results are tested. 

Subject headings: galaxies: evolution — galaxies: ISM — ISM: kinematics and dynamics — MHD — 
stars: formation — turbulence 
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1. Introduction 

The large-scale star formation (SF) cycle in disk 
galaxies plays a key role in modeling galaxy for- 
mation and evolution. This cycle crucially de- 
pends on the dissipativc properties of the turbu- 
lent interstellar medium (ISM). In normal, non- 
interacting disk galaxies, SF is believed to be 
statistically stationary, and self-regulated 1 by a 
balance between the energy injection and dissi- 
pation in the turbulent ISM in the disk (e.g., 
Scalo & Struck-Marcell 1984; Vazquez & Scalo 
1989; Dopita 1990; Firmani & Tutukov 1992, 
1994; Dopita & Ryder 1994; Wang & Silk 1994; 
Firmani, Hernandez, & Gallagher 1996; Avila- 
Reese 1998; Struck & Smith 1999; Avila-Reese 
& Firmani 2000). A self-regulating SF mecha- 
nism has also been used in semi-analytical mod- 
els of galaxy formation (e.g., White & Frenk 1991; 
Kauffmann, White, & Guiderdoni 1993; Cole et al. 
1994; Somerville & Primack 1999; van den Bosch 
1999), but in this case it has been applied to a 
hypotethic intrahalo medium in hydrostatic equi- 
librium with its own gravitational potential plus 
that of the cosmological dark matter halo, which 
extends up to approximately 15-20 times the opti- 
cal radius of the galaxy. In these models the disk 
ISM is virtually ignored and the feedback from 
stars is assumed to efficiently reheat and drive 
back the cooled gas into the intrahalo medium, the 
reheated gas fraction being a strong function of the 
halo mass. Unfortunately, in most galaxy forma- 
tion and evolution models, the detailed thermo- 
hydrodynamics of the ISM has not been treated 
explicitly. Therefore, several assumptions and ap- 
proximations have had to be made about the SF 
feedback and the dissipative properties of the ISM. 

Stars — mainly massive OB stars and the SN 
explosions they produce — are sources of thermal 
and turbulent kinetic energy in the ISM. Results 
from thermo-hydrodynamical simulations in dif- 
ferent contexts have shown that the kinetic en- 
ergy injected by SNe and massive stars to the 
ISM deeply affects the dynamics of the gas in disk 
galaxies causing it to be highly turbulent (e.g., 

Throughout the paper, by self-regulation we only mean 
that the SF rate feeds back on itself, generally in a neg- 
ative fashion, maintaining a statistically stationary value 
over long time scales, but of course allowing for significant 
spatial and temporal fluctuations. 



Bania & Lyon 1980; Chiang & Prendergast 1985; 
Chiang & Bregman 1988; Navarro & White 1993; 
Mihos & Hernquist 1994,1996; Rosen & Bregman 
1995; Vazquez-Semadeni, Passot & Pouquet 1995; 
Friedly & Benz 1995; Passot, Vazquez-Semadeni 
& Pouquet 1995; Gerritsen 1997; Avillez 1999; Ko- 
rpi et al. 1999; see also the reviews by Scalo 1987 
and by Vazquez-Semadeni et al. 2000). Cloud for- 
mation and disruption, the structure and dynam- 
ics of the vertical gaseous disk (in particular its 
scale height), fountains and chimneys, and even 
a huge hot gas corona (the "intra- halo medium") 
in virial equilibrium (including the potential of 
a huge cosmological dark matter halo), might be 
some of the phenomena and processes related to 
the turbulence produced by stellar kinetic energy 
sources. A key ingredient in all these phenom- 
ena and processes is the kinetic energy dissipation 
rate in the turbulent ISM, since it determines the 
effectiveness of large-scale SF feedback on galaxy 
formation and evolution. 

Dissipation in turbulent incompressible fluids 
is mainly controlled by a turbulent cascade from 
large to small scales, since it is in the latter 
where kinetic energy is dissipated. Several ana- 
lytical estimates and numerical simulations have 
shown that incompressible MHD turbulence de- 
cays as t~ a with 0.7<a< 1.0 (Biskamp 1994; Hos- 
sain et al. 1995; Galtier, Politano, & Pouquet 
1997). In the case of the large-scale ISM, turbu- 
lence involves supersonic MHD compressible non- 
isothermal flows, which most likely are dominated 
by shocks. High-resolution 3D simulations have 
been used recently to investigate dissipation in 
compressible MHD isothermal flows with param- 
eters corresponding to Galactic molecular clouds 
(Mac Low et al. 1998; Stone, Ostriker & Gammie 
1998; Padoan & Nordlund 1999; Mac Low 1999a). 
For the decaying regime, Mac Low et al. (1998) 
and Stone et al. (1998) have found E^ to decay as 
t~ a , with a ~ 0.8 — 1.0. For the forced regime, 
Stone et al. (1998) and Mac Low (1999a) con- 
cluded that the characteristic turbulent dissipa- 
tion time td is of the order of or smaller than the 
crossing time for the driving (or forcing) length 
at the rms turbulent velocity u rms , even in the 
presence of strong magnetic fields. In those sim- 
ulations, turbulence has been driven in Fourier 
space, with a fixed kinetic energy injection rate, 
E™, generating random large-scale velocity fluc- 
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tuations with a constant u rms . As a result, kinetic 
energy is injected ubiquitously (i.e., everywhere 
in space), and typically at large spatial scales, al- 
though Mac Low (1999a) also considered interme- 
diate injection scales. 

The way in which is injected into the ISM 
certainly differs from the prescription used in the 
studies mentioned above. The scales at which ki- 
netic energy injection occurs (the regions directly 
heated or accelerated by the stellar activity, of 
sizes of several pc to a few hundreds of pc), arc 
small compared to the scales of interest in the 
large-scale ISM (up to a few kiloparsecs). 2 Be- 
sides, the sources are of short duration and located 
at discrete, generally isolated sites. Moreover, the 
global energy input rate is not expected to be con- 
stant in general. This situation implies that in the 
ISM both forced- and decaying-turbulence regimes 
should coexist, the latter being the regions located 
far from the energy injection sources. 

In this paper we use two-dimensional (2D) nu- 
merical MHD simulations of turbulent compress- 
ible fluids resembling the ISM of normal disk 
galaxies (Vazquez-Semadeni et al. 1995, 1996; 
Passot et al. 1995) to explore the ability of the 
ISM to dissipate the turbulent kinetic energy in- 
jected by stellar sources. In order to do this, we 
use an ISM-like turbulent kinetic energy injection 
mechanism, based on "stellar winds" applied ran- 
domly in the medium; thus, we study how the 
dissipative properties of the flow depend on the 
parameters of the injection sources. 

It should also be noted that the ISM at large 
scales is far from isothermal (see, e.g., Myers 
1978), to the extent that this non-isothermality 
has constituted the starting point for multi-phase 
models of the ISM (e.g., Field, Goldsmith & 
Habing 1969; Cox & Smith 1974; McKee & Os- 
triker 1977). Therefore, in this paper we consider 
turbulent flows in the presence of parameterized 
radiative cooling and background heating (mim- 
icking the heating from a background UV field 
and cosmic rays). 

In §2 we describe the numerical method, the im- 

2 Strictly speaking, since such forcing is confined to a limited 
range of (small) spatial scales, in Fourier space it contains 
components of all spatial frequencies. In reference to the 
stellar driving we consider, in this paper the "scale of the 
forcing" will always refer to its characteristic size in physi- 
cal space. 



plcmentation of the kinetic energy injection mech- 
anism, and the parameters and initial conditions 
of the simulations. In §3 we discuss the dissipation 
time scales for driven (forced) turbulence simula- 
tions using various energy injection regimes, find- 
ing that most of the injected energy is dissipated 
locally; we also present a simple model to calcu- 
late id as a function of observable quantities of 
the ISM. In §4 we explore the decay and prop- 
agation of the "residual" turbulent energy that 
"survives" the strong local dissipation, by con- 
sidering the simulations after turning off the stel- 
lar energy input sources. In §5 we discuss pos- 
sible caveats of our simulations, such as the di- 
mensionality and the low resolution, and present 
higher resolution runs which support our conclu- 
sions. Then §6.1 compares our results with pre- 
vious work, and §§6.2 and 6.3 discuss the impli- 
cations of our results for models of galaxy evolu- 
tion and formation, respectively. Finally, in §7 we 
present a summary and our conclusions. 

2. Numerical method and simulations 
2.1. The model 

We use a slight variation of the numerical model 
presented in a series of previous papers (Vazquez- 
Semadeni et al. 1995, 1996; Passot et al. 1995) 
and the reader is referred to those papers for de- 
tails; here we just sketch the method. The MHD 
equations, including the internal energy conserva- 
tion equation, are solved in two-dimensions (2D) 
in the presence of model terms for radiative cool- 
ing, background heating, and stellar winds. While 
in previous papers self-gravity, the Coriolis force 
and large-scale shear were included, here we do not 
consider them because of various reasons. First, 
self-gravity causes runaway gravitational collapse 
when the stellar energy injection is not applied at 
the densest regions, as we do here (see §2.2 below). 
Second, the Coriolis force modifies the dissipation 
rate at long times in decaying regimes, but we con- 
sider this to be an unrealistic situation (see §4). 
Third, the large-scale shear introduces a non-zero 
lower bound to the kinetic energy, which we feel is 
irrelevant for the purposes of the present paper. 

The numerical technique is pseudo-spectral 
with periodic boundary conditions, and we use 
a combination of regular second-order viscosity of 
the form /i[V 2 u + (1/3)V(V • u)], and hypervis- 
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cosity of the form j/V 8 u in the momentum equa- 
tion. The hyperviscosity operator is a very steep 
function of wave number k, and, upon adequate 
choice of the coefficient v, confines viscous ef- 
fects to a much smaller range of scales than the 
regular second-order (or lower-than-cigth-order in 
general) viscosity, allowing a more extended in- 
crtial (self-similar) range in the kinetic energy 
spectrum. Our simulations show clear power- 
law inertial ranges over more than one decade 
in wavenumbers (§6.1). It should be pointed out 
that finite-difference schemes necessarily introduce 
numerical dissipation which is typically second-or 
third-order. On the other hand, hyperviscosity 
introduces spurious oscillations in the vicinity of 
strong shocks (Passot & Pouquet 1988). In or- 
der to "filter" these out we add a small amount 
of second-order viscosity, which, however, can be 
kept much smaller than what would be necessary 
without the hyperviscosity. Note that the use of 
hyperviscosity should not affect the global dissipa- 
tion rate of the flow, as this rate adjusts itself to 
counterbalance the energy injection by the stars, 
independently of the specific form of the dissipa- 
tion operator. What changes from one operator 
to another at a given dissipation rate is the dissi- 
pation scale (e.g., Frisch 1995) and the saturated 
energy content (Stone et al. 1998). 

Additionally, the code uses second-order mass 
diffusion (with coefficient /j, p ) with the purpose of 
smoothing out shocks even further, as the spectral 
scheme we use cannot handle exceedingly steep 
gradients. We do not expect the mass diffusion 
to affect the dissipation rates in the code since its 
only purpose is to prevent the formation of ex- 
cessively small density structures. The dissipative 
coefficients are chosen in such a way that the dissi- 
pation operators are already strongly active a safe 
margin before the shocks become too steep to be 
resolved. 

In the I28 2 simulations we use v = 1CP 11 , /i = 
0.003 and /j, p = 0.06. We take zero electrical resis- 
tivity, and a thermal diffusivity rj = 7[/x+^(n/4) 6 ], 
where 7 is the ideal-gas ratio of specific heats, and 
n is the resolution. All coefficients are expressed 
in the code's non-dimensional units, and are, in 
the momentum and magnetic flux equations, es- 
sentially inverse Reynolds numbers. 

We have chosen physical units in such a way 
that the simulations resemble the ISM in the plane 



of the Galaxy near the solar neighborhood at the 
1 kpc scale. We adopt as normalization values 
a number density n = 1 cm~ 3 , a temperature 
To = 10 4 K, a velocity uq = 11.7 km s _1 , and 
a length Lq = 1 kpc. The velocity and length 
units imply a time which (divided by 2ir since in 
non-dimensional units, the length of the box is 
2ir) gives the code time unit, to — (L /uo)/2n = 
1.36 x 10 7 yr. The magnetic field is written as 
B = Boe x + b, where Bo = 1.5/itG represents the 
uniform azimuthal component of the field and b 
is a superposed fluctuating component of rms am- 
plitude = 5fiG. 

2.2. "Star" formation and kinetic energy 
injection 

In previous papers, an "OB star" (or, in gen- 
eral, a source) was turned on at grid point x when- 
ever the density exceeded a certain (arbitrary) 
threshold value and the velocity field had finite 
negative divergence. A "star" consisted of a lo- 
cal source of heat, which caused the nearby gas 
to expand, forming warm bubbles resembling HII 
regions. After some experimentation, in this pa- 
per we have chosen to use "winds" (local outward 
radial accelerations, see Vazquez-Semadeni et al. 
1996) rather than heating, and to place the stars 
randomly with a given probability, rather than at 
the density peaks. 

The reasons for our choice are that the heating 
scheme deposited unspecified amounts of kinetic 
energy into the flow, the exact amount depending 
on the local cooling rate, which in turn depends 
on the local density and temperature, rendering 
the measurement of the kinetic energy injection 
rate impossible. Moreover, the density-directed 
placement of the sources caused three phenomena 
which obstruct the study of how the results de- 
pend on the the initial parameters of the energy 
injection mechanism. First, the OB star formation 
rate (SFR) fluctuated strongly, with alternating 
periods of dormancy and of strong activity (see 
Vazquez-Semadeni et al. 1995). This complicates 
the determination of the characteristic energy in- 
jection time, Ek/E™, where E™ is the global en- 
ergy injection rate, when the latter goes to zero. 
Second, the density-directed SFR produces strong 
clustering of the sources, because generally more 
than one grid point satisfy the density criterion 
within a density peak (a "cloud"). This causes 
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their outputs to interfere destructively (since op- 
positely directed accelerations acting on the same 
grid point cancel out) in such a way that the to- 
tal energy input of a cluster is smaller than the 
sum of the inputs from its constituent stars taken 
separately. Besides, due to clustering, large "su- 
pershells" form; these supershells are energy in- 
put sources which cannot be easily related to the 
initial injection parameters. In fact, destructive 
interference between the sources and "supershell" 
formation still occur in the random-star placement 
scheme at the largest SFRs, but to a much lesser 
extent. Finally, the density-directed scheme in- 
volves the presence of an initial negative velocity 
divergence at the forcing sites, making it impos- 
sible to determine the final velocity attained by 
the flow around the sources, since the acceleration 
acts on a pre-existing velocity field of fluctuating 
magnitude. In summary, the random placement of 
stars we use here, although certainly less realistic 
than the density-directed placement, makes for a 
much more controlled energy input, and we adopt 
it in our forced simulations. 

We add a fixed number of sources at each code 
timestep by specifying the probability with which 
a source is placed in each grid point. However, 
since the timestep is variable in the code, the in- 
stantaneous number of sources n s in the simula- 
tion is not perfectly constant in time. This is why 
we prefer to use the SFR per unit of area, Eob, 
measured directly from the simulations, in such a 
way that at any one time n s = Eob A^ob, where 
A^ob = 6.8 Myr is a typical lifetime of OB stars. 
Once a source has been turned on at a given grid 
point x, it stays on for a time AioB- During this 
time the gas around the source receives an accel- 
eration a directed radially away from x, with a 
magnitude 

a(r) = ^ exp[ _ (r 2_ (7 2 )/2(T 2 ]j (1) 

where a max is the amplitude of the acceleration, r 
is the distance from the position of the star, and 
a is a free parameter of the order of a few pixels. 
At a resolution of 128 grid points per dimension, 
1 pixel w 7.8 pc. The acceleration a produces an 
evolving velocity profile v(r, t) around the source. 
We show an example of the acceleration and veloc- 
ity profiles at the end of a stellar lifetime in fig. 1 
for the fiducial simulation called run 30 (see Table 
1). It is easily seen that a(r) has its maximum at 



r = a, and that by r — 2.2a, it has decreased to 
0.2a max . We choose this as a characteristic radius 
of the forcing region, which we denote If. It is im- 
portant to note that If is not the final radius of the 
expanding shells that result, which are actually a 
response of the flow to the energy injection mech- 
anism, and expand to sizes quite larger than If. 
This is particularly noticeable when stars cluster, 
in which case the shells may reach sizes of a few 
hundred parsecs, becoming "supershells" . 

The characteristic power ek and energy ek in- 
jected to the flow by one source can be roughly 
estimated as 

e k = I Z g u. ad i X * { Z g )fUfaf^ = MfUfaf 

JSi 

e k = e k Ai B, (2) 

where u is the flow velocity vector, Sf = nlf is 
the area of the forcing region, (£ g )f is its initial 
average surface density, and Uf and at are char- 
acteristic values of the velocity and acceleration, 
respectively. Note that, since the simulations are 
two-dimensional, it is necessary to specify some ef- 
fective disk height Hf for the forcing regions in or- 
der to define their surface density and mass. In the 
Galaxy, the regions of highest young cluster den- 
sity oscillate above and below the plane by ~ 50 pc 
(Alfaro, Cabrera-Cano, & Delgado 1991). There- 
fore, we take Hf = 100 pc. Then, (X g )t = poHf 
with po = ttih no, and Mf — Sf(£ g )f; we use 
the simulation mean density n as representative 
of the mean density in the forcing regions because 
the sources are seeded randomly in the simulation. 

The product afUf may be defined as the char- 
acteristic power per unit of mass injected per forc- 
ing region, e k . We define Uf = u{a,t — AtoB) < 
u m ax(t = AtoB) and Of as the value of a(r) 
where the area contained under the acceleration 
profile is divided horizontally in two equal re- 
gions. This roughly happens at r = 2a, where 
a(r) — 0.45a max ; therefore, we take at — 0.45a max . 
These values were decided upon empirically, but 
we have checked that they apply reasonably well 
to all of our simulations. Note that Uf is not triv- 
ially related to Of , AioB and if, due to the shape of 
the acceleration profile and to the fixed time over 
which the acceleration is active. In particular, we 
have checked that the simple formula v 2 ~ 2afZf 
does not apply very accurately. 

Together, the three parameters at, If and Uf 
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completely define the local energy injection pro- 
cess. 3 However, of and Uf have meaning mostly 
within the context of our simulations only, while 
oc U[Ci{l[Hf (or, cquivalently, ek, since we take 
AioB = 6.8 Myr=const in all cases) is a more 
physically and observationally relevant quantity. 
Thus, in what follows, we will consider ek and If 
as the determining parameters of the local prop- 
erties of the injection process. Additionally, there 
is a third, global, parameter, namely Sob- In 
the following sections we explore how the dissi- 
pative properties of the turbulent flow depend on 
the forcing parameters: If, ek, and Sob- 

2.3. The simulations 

We have performed a sizeable number of simula- 
tions, summarized in Table 1. Most of the simula- 
tions have been done at a low resolution of 128 grid 
points per dimension, in order to achieve reason- 
able coverage of parameter space. The dependence 
of the dissipation rate on resolution is further dis- 
cussed in §5, where test runs at higher resolution 
(512 2 ) are presented as well, for which the results 
obtained from the low-resolution simulations con- 
tinue to hold. 

We adopt the simulation labeled run 30 as the 
"fiducial" case. The set of simulations was initially 
chosen in such a way as to have at least two dif- 
ferent values of each of the three parameters If, Uf, 
and n s , for (roughly) constant values of the others 
(see Table 1). However, as our study progressed, 
we realized that the most relevant parameters are 
not exactly those we considered initially and, as 
a result, this philosophy is not clearly reflected in 
the set of derived parameters If, ek, and £ob- In 
any case, the three derived input parameters span 
a reasonable range of values and, as will be de- 
scribed in §3, we have constructed a simple model 
which combines them, to which our simulations 
will be seen to conform to within < 20%. 

In the fiducial run, the parameters of the ki- 
netic energy injection were chosen to represent a 
compromise between values typical of expanding 
HII regions and SN remnants at late evolution- 
ary stages (unfortunately, our code cannot handle 

3 We wish to stress here that while the former two parameters 
are controlled by us, the latter is already a response of the 
flow to the first two, and has to be measured from the 
simulations, rather than being user-defined. In this sense, 
our model here is still semi-empirical. 



very strong shocks). For run 30, the power and 
the total kinetic energy injected by each source is 
el k = 1-35 x 10 35 erg s" 1 and e k = 2.84 x 10 49 erg, 
respectively (see Table 1). The estimated radius 
of a typical HII region before the SN explosion 
is ~ 50 — 60 pc; at this time the kinetic energy 
deposited by the region -expanding with a veloc- 
ity of ~ 10 km s _1 on average- into the ISM is 
- 1 - 2 x 10 49 erg (e.g., Elmegreen 1991). The 
cooled shell swept up by a SN (type II) ejecta 
has roughly the same amount of kinetic energy, 
although this energy is injected in a time period 
much smaller than in the case of the HII region 
stage. 

We start the simulations with uniform density 
and temperature and nearly zero velocity every- 
where (u rms ~ 0.1km s _1 ), and begin placing 
sources randomly as described in §2.2. After one 
source lifetime, the number of sources in the sim- 
ulation reaches an approximately stationary num- 
ber. It is not precisely constant because of the 
probability-based algorithm we use and the vari- 
able timestep. The sources generate expanding 
shells which interact in a complicated fashion, pro- 
ducing a network of filaments and distorted shells, 
in a statistically stationary regime. Figure 2 shows 
the density and x-component of the acceleration 
at a typical instant in the evolution of the fiducial 
simulation, run 30. 

We have also performed a simulation with the 
density-directed scheme of source placement (run 
6) in order to compare our results to this more re- 
alistic simulation. The density threshold used for 
turning on a source was 4po- The input param- 
eters for this run were If = 30 pc and Uf =~ 7 
km s _1 (estimated from the final expansion ve- 
locity of the shells). Finally, in order to study a 
strictly decaying regime, we performed a simula- 
tion without SF activity (run 12 in Table 1). This 
run used a late evolutionary stage of run 6 as its 
initial condition. 

3. Dissipation in driven turbulence 
3.1. The global dissipation rate 

For all the forced runs described in §2.3 we mea- 
sure several global quantities of the flow: the ki- 
netic energy content, E^, the kinetic energy in- 
jection rate, E™, the total rate of change of E^, 
dEk/dt, and the rms velocity dispersion, u rms . 
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Figure 3 shows the evolution of these quantities 
for the fiducial run except M rm s, which we do not 
plot because, as shown in Fig. 4, it scales very 

1 /2 

closely as E k ' (see also Mac Low 1999a). These 
global quantities are seen to be statistically sta- 
tionary, without exceedingly strong fluctuations in 
time. The injection rate E™ is seen to fluctuate 
moderately, typically some 50% around its mean 
value. This is not only because the number of 
sources varies slightly in time, but also because of 
the complex interaction between the acceleration 
from the sources and the pre-existing velocity field 
(see eq. [1]). Similar behaviour is observed in all 
forced runs with random source placement. How- 
ever, a much more irregular behavior is present 
(not shown) in the run with density-directed SF, in 
which the SFR alternates between zero and strong 
bursts. Table 2 gives the average of these quanti- 
ties between t = 40.8 and 272 Myr for all runs. 

From the measured values of E k and E k n , we 
can define the energy injection timescale as 

t in = E k /E™. (3) 

However, the actual quantity of interest is id, the 
kinetic energy dissipation timescale, defined as 

id = E k /E^ s , (4) 

where E k lss is the energy dissipation rate. In a per- 
fectly stationary regime, E k and E^ lss coincide, 
as in the well-known Kolmogorov (1941) theory in 
the incompressible case, and in the simulations of 
Stone et al. (1998) and Mac Low (1999a) in the 
compressible case. However, in the more realistic 
case of a regime stationary only in an averaged 
sense, the two rates need not coincide. We thus 
estimate E k lss from the equation: 

where dE^/dt can be measured directly in the sim- 
ulation from the evolution of E k . Afterwards, id 
can be computed from eq. (4). The evolution of 
E k ms and i d is also shown in fig. 3. The latter 
quantity is shown averaged over periods of 20 suc- 
cessive outputs of the code, because of its high 
sensitivity to fluctuations in E™. The values of 
£^ iss and t d averaged between t — 40.8 and 272 
Myr are given also in Table 2. 



It should also be noted that eq. (5) is actu- 
ally incomplete, and should include transfer terms 
from (to) kinetic to (from) other energy forms, 
such as magnetic or thermal. Thus, such transfer 
terms are effectively included in £^ lss . In particu- 
lar, this includes kinetic energy losses due to first 
converting it to those other forms and then dis- 
sipating these by resistivity or radiative cooling. 
On the other hand, this implies that E^ iss may 
also act as a source of kinetic energy, in which case 
it may become negative. However, it is seen that 
^diss < q Qn iy during a brief period at the begin- 
ning of the simulation, which, together with the 
similarity between E™ and E^ lss discussed below, 
reassure us that the error incurred by omitting the 
transfer terms in eq. (5) is not too severe. 

Remarkably, fig. 3 shows that E k and E* iss 
are always very close to each other. Since the en- 
ergy injection occurs at very localized and scat- 
tered sites, while the dissipation rate is computed 
over the whole flow, their highly similar evolutions 
suggest that most of the dissipation must occur 
at or near the injection regions, and dominates 
the global dissipation rate 4 . Only in this case 
can the global dissipation rate track the locally- 
originated injection so closely. In turn, this implies 
that only a small fraction of the energy injected at 
the sources "escapes" to more distant regions. We 
refer to this as the "residual" turbulence, which 
coexists in the flow with the forced regions. We 
discuss the propagation of turbulent kinetic en- 
ergy from the injection to the quiescent sites in 
§4. As a consequence, in the remainder of this pa- 
per we refer to the injection and the dissipation 
times indistinctly. 

Analyzing all the forced runs, we find that id 
(ii n ) does not significantly depend on two of the 
energy injection parameters, namely the energy in- 
put rate of the source, e^ (or just its total energy 
input e k ) and the source formation rate, Sob (see 

4 A posteriori, this is actually not surprising; in stationary 
forced incompressible regimes, it is well known that the dis- 
sipation rate is proportional to (see, e.g., Frisch 1995; 
see also Stone et al. 1998; Mac Low 1999a for the com- 
pressible case). Since the kinetic energy density is greatest 
at the injection sites, most of the dissipation must occur 
there as well. Thus, the dissipation rate is actually a func- 
tion of position (through the local kinetic energy density) 
and may fluctuate strongly from one position to another, 
as is also the case in incompressible turbulence (see, e.g., 
Frisch 1995). 
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§2.2 for definitions). The relevant forcing param- 
eter seems to be the forcing scale, If. We obtain 
the following approximate empirical correlation of 
id on the input parameters of the forcing: 



;0.60 

'f 



;0.84 

'f 



(e k Sqb) 012 (e k Sob) 012 



(6) 



3.2. Relating t- ln to measurable properties 
of the ISM 

In order to estimate theoretically £d(~ tin) as 
a function of frequently used astronomical quanti- 
ties of the ISM, we can use the definition of ti n : 



/ A £ g u 2 dA' 



(7) 



Jk / A S g a-udA'' 
The integral in the numerator can be estimated as 



H g u 2 dA' « <S g ) u 2 ms A, 



(8) 



where (S g ) is the average gas surface density and 
A is the total area where the flow is contained. 
Similarly to what was done in eq. (2), the integral 
in the denominator can be estimated as follows: 

/ S g a • udA' = / S g a • udA' rts (S g ) f a f Uf Af 
J A J A t 

= (£g>f«Mf, (9) 

where the first equality follows because a = out- 
side of the forcing regions, whose total area is Af. 
Therefore, from eqs. (7), (8), and (9) we obtain 



4/ 



(10) 



where / = Af/A is the filling factor of forcing re- 
gions. Naively, one could take / = irl 2 n s /A. How- 
ever, one has to take into account the fact that 
some of the forcing regions overlap; the probabil- 
ity of this occurring is larger at larger if. Thus, 
f <irl 2 n s /A in general, and in practice we prefer 
to measure / directly from the simulation as the 
fraction of the simulation area in which the stel- 
lar acceleration is larger than 20% of a max (see eq. 
[1]), corresponding to the area within a distance If 
from the centers of the forcing. 

Table 2 gives the values of u rms and / measured 
in all the forced runs (e k is an input parameter 



given in Table 1), as well as the values of U n re- 
sulting from them. These "theoretical" injection 
times are seen to be reasonably close to those mea- 
sured directly in the simulations according to eq. 
(3). For most of the runs, the agreement is better 
than 20%, with the largest deviations being ~ 45% 
and 30%. 

A posteriori, it is not surprising that eq. (10) 
gives estimates for t[ n consistent with the values 
obtained using eq. (3), since, after all, we have 
done nothing more than estimating the integrals 
intervening in eq. (7) in terms of averaged quan- 
tities. However, the relevance of this result is 
twofold. First, it shows that reasonable estimates 
can be obtained using global averaged quantities 
even when the dissipation is extremely localized, 
and second, it provides a means of relating the dis- 
sipation time to frequently used quantities of the 
ISM. Indeed, using the expression f <itl 2 n s / A for 
the filling factor in eq. (9) we obtain 



S g a • ugL4'^ (E g )fe k 7Wf n s = 
Mfe k MoB^OB,A = e k S B-4, 



(11) 



where Mf — (S g )f7rZ 2 , e k = Mfi^AtoB, and Sob 
has been defined in §2.2. Using eqs. (8) and (11) 
in eq. (7), one obtains 



iin> 



<S> 



2 

rms 



e k S Q B 



(12) 



which is our desired connection to measurable 
quantities of the ISM. For example, taking w rms = 
8 km s -1 , S b = 3 x 10~ 5 yr~ 4 kpc~ 2 (Tammann, 
Loffler & Schroeder 1994), e k = 3.5 x 10 49 erg, and 
the gas surface density typical for the solar vicin- 
ity, S g w 12M pc- 2 , we find that t in > 1.5 x 10 7 
yr. Note that the value of e k we used reflects a 
contribution of roughly 1.5 x 10 49 ergs per HII re- 
gion plus ~ 2 X 10 49 ergs of mechanical energy per 
SN event. The resulting value of t ln is in rough 
agreement with our fiducial simulation. 

It is worthwhile to note that t- m depends on 
S g and Sob, and since both of these vary with 
Galactocentric distance, one expects t ln (and con- 
sequently the dissipation time) to exhibit such a 
dependence. This question will be discussed in 
further detail in section §6.2. 

One more point should be remarked regarding 
eq. (10) (or, equivalently, relation [12]), namely 
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that the numerical simulations clearly show that 
the dissipation time increases with the size of the 
forcing regions (see eq. [6]), even though relation 
(12) has no explicit dependence with If. This is 
equivalent to asking the question: how does the 
empiric expression for i; n (cq. [6]) connect with 
relation (12), which was derived directly from the 
definition of t ln ? These two relations imply that 
u rms must scale approximately as If ■ 42 (ekSoB) ' 44 - 
Recalling that in relation (12) the inequality ac- 
counts for overlapping of the regions, one expects 
that an extra factor depending on If must be in- 
cluded in the RHS of relation (12) in order for it to 
become an equality. In the simulations we empiri- 
cally find that roughly u rms cx Z f - 40 (e k S O B) - 42 . A 
full MHD theory of the problem should account 
for this dependence. 

3.3. A more realistic simulation 

In the runs presented above, we have not in- 
cluded self-gravity nor the Coriolis force, and the 
energy input sources were placed randomly with a 
fixed probability, in order to have good control of 
the forcing parameters, and to uncover the depen- 
dences presented in §§ 3.1 and 3.2, at the expense 
of some realism. Somewhat more realistic sim- 
ulations should include those additional agents, 
and involve a density-directed placement of the 
sources. In particular, the latter scheme causes 
the energy injection rate to fluctuate strongly and 
the input sources to be highly clustered; large su- 
pershells form and significant cancellation of the 
energy input from neighboring sources occurs (see 
§2.2). It is thus interesting to see if our results still 
hold in this case. 

To this effect, we have also performed a simula- 
tion with these ingredients (run 6, see §2.3). In fig. 
5 we present images of density field of this simula- 
tion at times t = 14.3 Myr (left) and t = 144.3 Myr 
(right). At the earlier time, small shells arc seen, 
whose sizes are comparable to the sphere of influ- 
ence of the "winds" themselves (of radius If = 30 
pc). At the later time, the shells are seen to be 
larger, but their larger sizes are not only a con- 
sequence of the inertial motions induced by the 
forcing, but also of induced-SF events in the shells 
themselves. Even though for this run we cannot 
estimate Uf precisely, and consequently ek, we can 
measure the total injection rate and, from simple 
visual inspection, we see that the average radius 



of the shells in the simulation is at least twice as 
large as 2lf. Moreover, on average we measure 
Ek w 0.18 (in code unities), u lms w 7 km s _1 , and 
td w 18 Myr. We thus see that this run is in a 
very similar regime as runs 30 and 35, and in par- 
ticular the dissipation time is within 15% of that 
measured in these runs, and is 30% larger than 
the lower-bound estimate for the solar neighbor- 
hood given in §3.2. Thus, this run reassures us 
that the results obtained with the random place- 
ment of the stellar sources are applicable in more 
realistic situations. 

4. Dissipation in decaying turbulence 

In order to study the decay of turbulence far 
from the injection sites, we now consider a simula- 
tion in a decaying (i.e., without energy injection) 
regime. To this end, we consider a continuation 
of run 6, named run 12, but with the SF turned 
off. In fig. 6 we present images of the density field 
from this run at two times, one immediately af- 
ter turning off the SF (t = 0) and the other at a 
much more advanced stage (t — 180.7 Myr). The 
former time shows a large number of expanding 
shells, while the latter contains no shells anymore. 
The evolution of the kinetic energy for this simu- 
lation is shown in fig. 7. We find that E k decays 
roughly as 

E k (t)=E {l + ^y a (13) 

with a « 0.8, in reasonably good agreement with 
previous studies of decaying MHD turbulence for 
isothermal flows (Mac Low ct al. 1998, Stone et 
al. 1998), and t is the code time unit. A charac- 
teristic decay time can be defined in this case as 
the time E k to decay by a factor of two at the be- 
ginning of the simulation. For a w 0.8, this occurs 
at t = 1.38 x t = 18 Myr. 

We can explore the decay of turbulent motions 
somewhat further. This analysis is greatly simpli- 
fied in the case of long times (» to), which allows 
us to drop the term 1 inside the parenthesis in eq. 
(13), writing 

E k (t)^E ^y a (14) 

This limit is reasonable, since we are interested in 
the behavior of the turbulence far from the injec- 
tion sites. Then, assuming that the rms velocity 
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u rms scales roughly as the square root of the total 
kinetic energy (see Fig. 4), from cq. (13) we can 
write: 

/j \ -a/2 
U Ims (t) W u i[^- ) . ( 15 ) 

where U\ is the rms velocity dispersion at t = to 
(?« 7 km s _1 ). Similarly to the forced case, here we 
may formally define the dissipation time scale as 
the ratio between and \Ek\- Thus, we obtain: 

\E k \ a\t J 

Thus, the characteristic dissipation time increases 
linearly with time in the decaying case, implying 
that the instantaneous dissipation time of the flow 
depends on the net energy content of the flow, be- 
ing smaller for larger energy contents. 

Identifying now the decaying case with the re- 
gions far from the injection sites in a forced case 
with small-scale, spatially intermittent forcing, the 
temporal decay of the energy can be translated 
into a decay with distance. We can estimate the 
distance "reached" by the turbulent motions in a 
time t as 

£(t) ~u rms (t)t «4(iy~ Q/2 , (17) 

where £o = Wito ~ 10 2 pc is the typical length 
traveled by a fluid parcel at the characteristic 
time and velocity at the beginning of the decaying 
regime. This length does not exceed the size scale 
of the expanding shells produced by the forcing. 
Solving for t in eq. (17) and substituting in eqs. 
(14) and (15), we obtain the dependence of E k and 
ii ms with £: 

E^l)^E {£)~^ (18) 

u rms (f)~ui(-) (19) 

For a w 1 this implies that E^ and u Tms decay 
with distance as (£/£o)~ 2 and (^/£o) _1 , respec- 
tively. Thus, "residual" turbulent motions with 
an rms velocity of roughly 6-10 km s _1 , as ob- 
served in our and other disk galaxies, would decay 
to roughly 1/10 of that value on distances ~ 0.6 — 1 
kpc, in the absence of energy sources. Note how- 
ever that, since the decay law (14) is a power law 



(rather than, say, an exponential), the time re- 
quired for Ek to decrease by consecutive constant 
factors increases as the the kinetic energy content 
of the flow decreases (cf. eq. [16]). This applies 
also to the decay with distance, and thus the dis- 
tance required for the rms velocity to decrease by 
a factor of 10 depends on the initial energy con- 
tent of the medium. This is also why energy is 
dissipated so efficiently near the injection sites. 

5. Effects of resolution and dimensionality 

For economy, and to maximize our ability to 
explore parameter space, in this paper we have 
limited ourselves to low-resolution, 2D numerical 
simulations. However, even though in this paper 
we are not attempting to give detailed quantitative 
results, but only first-order estimates, it is impor- 
tant to test that even these can indeed be attained 
with the simulations we have employed. 

The two-dimensionality poses an apparently se- 
rious problem, since it is well known that, in the 
incompressible case, the energy dissipation rate 
has a fundamentally different behavior in 2D with 
respect to 3D. While in 3D this rate is believed to 
remain finite even in the limit of vanishing viscos- 
ity, in 2D it does tend to zero in this limit (see, 
e.g., Lesieur 1990, sec. IX. 3. 5). Moreover, in 2D 
the kinetic energy cascades in an "inverse" way, 
from small to large scales (always in the incom- 
pressible case). All of this questions the validity 
of using 2D simulations for estimates of the dis- 
sipation rate and, in particular, implies that con- 
vergence tests may be meaningless in 2D since the 
kinetic energy dissipation rate does not converge 
as the resolution is increased and the dissipation 
coefficients are decreased. 

Fortunately, this problem may not be present in 
the highly compressible case, in which shocks dom- 
inate the dissipation rate, since it is well known 
that shocks cause direct transfer from all scales to 
the dissipation scales (Kadomtsev & Petviashvili 
1973). This process is independent of the dimen- 
sionality of the flow since the shocks are essen- 
tially one-dimensional structures. Moreover, the 
dissipative nature of the shocks suggests that the 
problem of vanishing dissipation in the limit of 
vanishing viscosity in 2D is eliminated in their 
presence. This suggestion is supported by the fact 
that in our forced simulations the dissipation rate 
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follows the injection rate closely. In summary, we 
are confident that in the highly compressible 2D 
case, the dissipation rate can still be meaningfully 
measured, and that convergence tests can be per- 
formed. 

To this effect, we have produced a higher- 
resolution simulation in the decaying regime, at 
512 2 grid points, still in 2D, labeled run 17. In this 
run, we have reduced the hyperviscosity coefficient 
by a factor of 4 7 = 16384, and the second-order 
viscosities and diffusivities by a factor of 4 1 = 4, 
in an attempt to cause the dissipative and diffu- 
sive scales to be a factor of 4 smaller than in the 
128 2 runs, i.e., roughly the same size in grid points. 
The factors chosen take into account the order of 
the hyperviscous and second-order viscous opera- 
tors (8 and 2, respectively), and the fact that the 
mean Fourier amplitude of the velocity scales as 
fc _1 if the kinetic energy spectrum E(k) scales as 
fc~ 2 , as expected in a flow dominated by shocks 
(see, e.g., the review by Vazquez- Semadeni et al. 
2000 and references therein). As is the case with 
the decaying run at 128 2 resolution (run 12), run 
17 is again a restart of a forced simulation which 
used the old scheme of placing the stars at the 
highest density sites. 

Figure 8 shows the evolution of the kinetic en- 
ergy Ek for run 17. It is seen that over the early 
epochs of the simulation (£ = log(l + t/to) ~ 0.4), 
the energy varies nearly as a power law with time, 
with exponent ~ —0.9, in rough agreement with 
the —0.8 exponent observed for run 12, and also 
with the exponents reported by Mac Low et al. 
(1998) and Stone et al. (1999). However, at later 
times, a departure from this behavior is observed, 
with the decay rate of E^ monotonically decreas- 
ing for 0.4<£<0.9, indicating that at the higher 
resolution of run 17 the decay rate indeed becomes 
much smaller after the shocks have subsided. To 
support this interpretation, in the same figure we 
also show the evolution of the most negative value 
of the divergence of the velocity field (shown in log 
of its absolute value). It is seen that the departure 
from a power-law decay rate coincides with a drop 
of the divergence by nearly a factor of 3 at £ ~ 0.5. 
Moreover, a sharp drop in Ey at £ ~ 0.95 coincides 
with a large peak in the divergence, after which 
the latter decreases even further and returns 
to a slow decay rate (£ ~ 1.1). In summary, Fig. 8 
supports the view that in 2D the dissipation rate 



follows roughly the same decay rate as 3D simula- 
tions (Mac Low et al. 1998; Stone et al. 1999) as 
long as the shocks dominate the dissipation, giv- 
ing us confidence that our 2D simulations provide 
realistic measurements of the dissipation rate in 
the forced regime, which contains large numbers 
of shocks, and in the early epochs of the decaying 
regimes. 

Concerning the forced potential problem 

is that for some of the low-resolution runs, the in- 
jection scale overlaps with the dissipation scale (J. 
Brasseur 2000, private communication), and thus 
there is the risk that this energy never reaches the 
larger scales of the flow. We believe this is not 
really a problem, because, as mentioned above, 
shocks "short-circuit" the energy directly from the 
energy-containing scales to the dissipations scales, 
so the cascade through intermediate scales is elim- 
inated anyway. In physical space, what happens 
is that, even if the forcing has characteristic scales 
larger than the dissipation scale, its effect is to 
induce shocks into the medium, which are small- 
scale structures. Indeed, runs 40 and 41, which 
have a value of If twice as large as that of run 30, 
still agree with our model prediction within ~ 45% 
and 30%, respectively. 

Nevertheless, it is important to verify that the 
agreement between our semi-empirical model of 
§3.2 and the low resolution simulations is main- 
tained in higher-resolution runs. To this effect, we 
have performed a 512 2 forced run, labeled run 44 
(Table 1) in which the physical size of the box (1 
kpc) is maintained, but the smallest resolved scale 
is now 1.95 pc. This run uses the same viscosities 
and difussivities as run 17, and its parameters are 
chosen so that it is a rescaling of run 30 (the fidu- 
cial run): If is 4 times larger in pixels, so it keeps 
the same physical size in parsecs. However, the 
dissipation scale is four times smaller in physical 
size, becoming safely smaller than the injection 
scale. All other parameters are the same as in 
run 30. From Table 2, it can be seen that run 44 
conforms nicely to our semi-empirical model, its 
predicted and measured injection times agreeing 
within ~ 15%. 
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6. Discussion 

6.1. Comparison with previous work 

The main difference between our MHD simu- 
lations of supersonic compressible flows and pre- 
vious ones aimed at exploring turbulence dissipa- 
tion resides in the method we use to drive the tur- 
bulence. As pointed out in §2.2, we have tried 
to represent the way in which kinetic energy is 
injected into the ISM of disk galaxies by stellar 
sources. Moreover, since our interest in this pa- 
per focuses in the large-scale ISM, we have taken 
parameters of the flow corresponding to the solar 
neighborhood, as done in previous papers model- 
ing the ISM at large (Chiang & Bregman 1988; 
Rosen & Bregman 1995; Vazquez-Semadeni et al. 
1995, 1996; Passot et al. 1995), and have not as- 
sumed the flow to be isothermal. All of this con- 
trasts with recent studies of turbulence dissipation 
in molecular clouds (Mac Low et al. 1998; Stone 
et al. 1998; Padoan & Nordlund 1999; Mac Low 
1999a), which have considered isothermal flows 
with large-scale, ubiquitous forcing. On the other 
hand, those works have been based on 3D simu- 
lations, while here we have used 2D simulations. 
However, large differences on the dissipation rates 
and timescales between 2D and 3D simulations are 
not expected (see e.g., Stone et al. 1998 and §5 
above). Specifically, the kinetic energy dissipation 
rates in 3D are expected to be slightly higher than 
in 2D, particularly in the regime of decaying tur- 
bulence, because at higher dimensionality, more 
shocks are produced due to the larger number of 
degrees of freedom (e.g., Mac Low et al. 1998). 
For example, Stone et al. (1998) reported that 
the energy decay times found for 2 + 1 / 2D MHD 
simulations are a factor of 1.50 — 1.75 times larger 
than those for 3D. 

At the results level, in spite of the very different 
physical regimes we consider, we agree with previ- 
ous works in that compressible MHD flows dissi- 
pate their turbulent kinetic energy very efficiently. 
More specifically, for their molecular cloud simu- 
lations, Stone et al. (1998) and Mac Low (1999a) 
find that the dissipation time is of the order to the 
flow crossing time at the energy-containing scale 
Ae, i.e. ^ « AE/u rms . However, in their work, 
this scale is very close to the injection scale be- 
cause of the large-scale and ubiquitous nature of 
the injection. In our case, on the other hand, Ae is 



in general quite larger than the driving scale 2Zf , 
because of the expansion of the shells, which is 
a way of transferring the energy to larger scales, 
and possibly also because energy that is not dis- 
sipated in shocks may escape towards the largest 
scales through an inverse cascade due to the two- 
dimensionality (to be investigated elsewhere). In- 
deed, fig. 9a shows the energy spectrum for run 
30 at times t = 1.36, 6.8 and 136 Myr. For ref- 
erence, fig. 9b shows the x-component of the ve- 
locity field at t = 6.8 Myr. It can be seen that 
at t = 1.36 Myr, the spectrum contains only one 
peak at log A: ~ 1.4 (the bulk of the spectrum at 
lower k at this time corresponds to the very mild 
initial conditions with u rms ~ 0.1 km s _1 ). A scale 
(in pixels) corresponding to wavenumber k can be 
defined as A = 128/fc. For logfc = 1.4, A 5 pix- 
els, and indeed at this time there is only one ex- 
panding shell in the simulation with "radius" ~ 4 
pixels, measured as the distance from the center 
to the velocity maximum of the shell (this dis- 
tance has been measured from cuts through the 
x-component of the velocity passing through the 
shell center, similar to those in fig. 1, not shown). 

At t = 6.8 Myr, this peak has moved to 
logfc ~ 1.22 (A w 7.7 pixels) due to the expan- 
sion of the shell, whose radius has grown to <~ 6 
pixels. New peaks at larger k have appeared cor- 
responding to the formation of new sources. The 
maximum of the spectrum (roughly, the "energy- 
containing scale") at t — 6.8 Myr, at logfc ~ 0.8, 
corresponds to the diameter of the largest shell at 
that time. Energy at even larger scales must cor- 
respond to the size of the region spanned by the 
various sources at that time, since the shells them- 
selves have not grown further than log fc = 0.8 yet. 
By t = 136 Myr, the spectrum is very close to a 
power law in the range 0.6 < logfc< 1.5, suggesting 
the existence of fully developed turbulence at that 
time. 

Since the results above suggest a clear separa- 
tion between the injection and energy-containing 
scales in our simulations, in order to compare our 
results to the works of Stone et al. and of Mac 
Low, we can quantitatively estimate the energy- 
containing scale in our simulations as the scale 
corresponding to the centroid of the energy dis- 
tribution, given by (see, e.g., Novikov 1978; San- 
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tangelo et al. 1989; McWilliams 1990) 

where Uk is the Fourier amplitude of the mode 
with (vector) wavenumber k. Table 2 gives the 
values of Ae for all runs at t = 136 Myr, a time 
at which the spectrum has become essentially sta- 
tionary, and the value of the large-scale crossing 
time t CT = As/urnis, calculated with the rms ve- 
locity at that same time. It is seen that the cross- 
ing times do not correlate well with the measured 
dissipation times, and are in general quite larger. 
We speculate that this is due to the fact that i cr 
measures essentially the dissipation time at the 
largest scales, and fails to capture the local dis- 
sipation near the sources. In particular, the lo- 
cal and large-scale dissipations contribute in dif- 
ferent proportions to the total rate depending on 
the size, number and specific energy injection rate 
of the sources, and thus the crossing time is not 
a good estimator of the true injection and dissi- 
pation times in the case of non-ubiquitous, small- 
scale forcing. 

In conclusion, we do not expect that a simple 
law of the type id = ^e/U, with U either the rms 
of the characteristic forcing velocity, can hold in 
the case with non-ubiquitous forcing: another pa- 
rameter, the filling factor of the forcing regions, is 
crucial, as it determines the total amount of en- 
ergy injected by sources of given local parameters, 
the proportion in which Uf contributes to u rms and 
the amount of overlapping between the regions. 

6.2. Implications for models of galaxy evo- 
lution 

From the point of view of galactic evolution, 
the global ISM processes, which take place at large 
space and time scales, are the most relevant. As 
mentioned in the Introduction, several galactic 
evolution models are based on the idea that the 
SF rate is controlled by a disk vertical balance in 
the ISM between the rates of energy injection due 
to SF and energy dissipation. Since turbulent ki- 
netic energy is believed to be one of the dominant 
sources of pressure in the ISM, the relevant dissi- 
pation rate for these models is that of turbulence, 
motivating our study. We have shown that is 
dissipated very efficiently near the input sources 



(locally). This implies that very large "active" 
turbulent zones can exist in galaxies only if the 
forcing regions are themselves large, as a conse- 
quence of either very powerful sources, or strong 
clustering and self-propagating SF. This might be 
the rule for starburst galaxies, but generally the 
exception in normal disk galaxies. Besides, as re- 
sults from our simulations have shown (§3), the 
dissipation time id indeed increases as the source 
size and filling factor increase. 

In the energy balance in the vertical direction 
assumed in the galaxy evolution models mentioned 
above (taking a given energy input rate, given by 
the SF rate, and estimating the dissipation rate 
as cx uj? ms /id), the equilibrium gas velocity disper- 
sion is determined by id- Hence, assuming ver- 
tical (one-zone) hydrostatic equilibrium, the disk 
HI thickness will depend basically on id- Accord- 
ing to the results of these galaxy evolution models 
(e.g., Firmani et al. 1996; Avila- Reese & Firmani 
2000; Firmani & Avila-Reese 2000), a realistic HI 
gas thickness at the solar neighborhood is obtained 
with values of id roughly twice as large as the val- 
ues found here (all the other relevant quantities as 
gas surface density, circular velocity and star for- 
mation rate, are self-consistently calculated and 
agree well with observations). From the hydro- 
static equilibrium condition, indeed it is easy to 
see that a turbulent gas layer supported by a veloc- 
ity dispersion « 9 km s _1 would have a thickness 
roughly 3 times smaller than observed, suggesting 
that probably other kinds of pressure, such as the 
magnetic field and cosmic-rays, are also important 
for supporting the galactic gas layer. 

On the other hand, it is important to have in 
mind that the vertical galactic disk is actually 
stratified; as the gas density decreases, the sur- 
vival and propagation of turbulent motions could 
be easier, and therefore id would be longer. In 
fact, shells are well known to accelerate when they 
transit from a high- to a low-density medium. In 
this sense, the dissipation times calculated here for 
a fixed gas average density can be interpreted as a 
lower limit; hence, the disk thicknesses estimated 
with our turbulent dissipation times will be also a 
lower limit. In any case, according to the results 
of galaxy evolution models, with self-regulated SF, 
the dissipation time obtained here is not too far 
from the values necessary to reproduce several ob- 
servational data. Thus, our results do not disagree 
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with the idea that turbulent motions induced by 
SF feedback contribute significantly to support the 
vertical HI disk (e.g., Lockman & Gehman 1991; 
Firmani & Tutukov 1992, 1994; Ferrara 1993; Do- 
pita & Ryder 1994), although some extra support 
may be necessary to keep the HI and gas ionized 
layers up at the observed heights. 

Finally, it is important to remark that the 
turbulent dissipation time may change along the 
galactic disk since the gas densities and dynami- 
cal conditions change with radius. According to 
Firmani et al. (1996), t^ oc R when the rotation 
curve is flat. In §3.2 we obtained a theoretical 
expression for t^ (actually ti n ; eq. [12]) which in 
particular implies a dependence on galactic radius 
through (Eg) and Sob- The OB-star formation 
rate is directly proportional to the total formation 
rate and, in most disk galaxies, oc £™, with 
n « 1.4 — 2.0 (e.g., Kennicutt 1998 and references 
therein). Therefore, from eq. (12) one obtains that 
tin(r) oc (S g )" _1 (r). The gas surface density typ- 
ically decreases with radius. For example, in the 
Galaxy, from 4 to 16 kpc, S g decreases roughly 
as r _1 (Dame 1993). Hence, the dissipation time 
should decay as the Galactocentric radius R to 
some power, probably less than unity. 

6.3. Implications for models of galaxy for- 
mation 

Several approaches of galaxy formation within 
the context of the hierarchical CDM-based sce- 
nario — in particular the so-called semi-analytical 
models (e.g., Kauffmann ct al. 1993; Cole et al. 
1994; Somervillc & Primack 1999; van den Bosch 
1999) — require the feedback from SF to be able to 
reheat the cooled gas up to the virial temperature 
of the whole dark matter halo. This requirement 
implies that a significant fraction of the released 
SN energy (probably more than 10%) remains as 
kinetic energy and that this energy is not dissi- 
pated within the disk. 

The physical picture invoked in the semi- 
analytical models implies a self-regulated SF not 
at the level of the disk ISM, but at the level of the 
hypothetical intra-halo medium. The SNe formed 
in the disk are assumed to reheat the cold disk 
gas up to the virial temperature of the cosmo- 
logical dark matter halo, driving it back into the 
intra-halo medium, and occasionally expelling it 
completely from the system (White & Frenk 1991). 



Thus, a crucial question for models of galaxy for- 
mation is whether the energy released by SNe and 
stars in the disk is able not only to maintain the 
warm and hot phases and the stirring of the disk 
ISM, but also to maintain a huge hot gas corona 
in quasi-hydrostatic equilibrium with the cosmo- 
logical dark halo (the intra-halo medium); the es- 
timated sizes of these dark halos are several tens 
and hundreds of kiloparsecs for dwarf and giant 
galaxies, respectively. Note that this intra-halo 
medium should not be confused with the diffuse 
ionized gas and high- velocity-dispersion HI gas at 
scale heights of one or a few kiloparsecs above the 
disk plane, also frequently referred to as a "halo" 
in the ISM community (see e.g., Reynolds 1997; 
Kalberla & Kcrp 1998; Mac Low 1999b). 

The question is then whether the intra-halo gas 
can be dynamically heated by the turbulent en- 
ergy input due to stellar winds produced by SNe 
and OB stars. Indeed, numerical simulations show 
that up to 5 — 10% of the energy produced by SNe 
can be transformed into kinetic energy in the ex- 
treme case of multiple SN explosions forming su- 
pershells (e.g., Silich et al. 1996; Korpi et al. 1999); 
for isolated SNe, this fraction is lower than ~ 1%. 

However, according to the results obtained 
here, most of kinetic energy in these "active" tur- 
bulent regions is dissipated locally and the w rms 
of the "residual" turbulent motions quickly de- 
cays with distance (roughly u lms oc i^ 1 ). In order 
for the turbulent gas to be driven back into the 
intra-halo medium, the typical sizes of most of 
the forcing shells should exceed the gaseous disk 
height (supershclls). This could be a common sit- 
uation in dwarf starburst galaxies, but it is not 
the rule in normal disk galaxies. In fact, even 
in dwarf starburst galaxies, gas ejection from the 
disk seems to be not very efficient, according to 
recent hydrodynamical simulations by Mac Low 
& Ferrara (1999). Besides, numerical simulations 
have shown that the disk magnetic field is an ef- 
ficient shield able to confine most superbubbles 
preventing immediate blowout (e.g., Slavin & Cox 
1992; Franco et al. 1995; Tomisaka 1992, 1998). 
Therefore, we conclude that it is very unlikely 
that the turbulent kinetic energy injected by SNe 
and massive OB stars into the disk ISM can sur- 
vive, propagate and drive back the gas into the 
intra-halo medium as the semi-analytic models of 
galaxy formation require. 
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We should note that our results are for densi- 
ties typical of the ISM in the disk plane at the 
solar neighborhood. One may think that the disk 
stratification might reduce the dissipation rate so 
much that the turbulent kinetic energy could actu- 
ally reach the intra-halo medium. However, since 
the HI disk scale height is larger than the typi- 
cal dissipation lengths we have found, we believe 
any escaping energy will be too weak to provide 
support to the intra-halo medium. Simulations in 
a stratified medium will be presented elsewhere, 
with a code that does not require periodic bound- 
ary conditions. 

An important question not treated in this paper 
is where the turbulent dissipated energy goes. It 
is well known that shocks efficiently transform ki- 
netic energy into thermal energy. Thus, the heated 
gas at relatively large heights could be a source of 
high-energy photons able to ionize the low den- 
sity extraplanar medium or even the hypothetical 
intra-halo medium. For example, Slavin (2000) 
has showed that hot gas in cooling SN remnants 
can be an important source of photoionization for 
the diffuse ionized gas. Dissipation in turbulent 
mixing layers has been also proposed as a mecha- 
nism to produce energetic radiation (Slavin, Shull, 
& Begelman 1993). Minter & Spangler (1997) 
have found that turbulent dissipation heating like 
that due to the ion-neutral collisional damping in 
the fluid-like turbulence model is able to provide 
a substantial contribution to the energy budget of 
the diffuse ionized gas. More studies on dissipa- 
tion of turbulence and the thermo-hydrodynamics 
of the extraplanar and intra-halo medium are cer- 
tainly necessary. 

7. Concluding remarks 
7.1. Summary 

We have used 2D simulations of MHD, com- 
pressible, self-gravitating flows in the presence of 
parameterized heating and cooling and wind-like 
stellar kinetic energy input in order to study tur- 
bulence dissipation in the ISM. The "winds" are 
described by the typical radius If of their region of 
influence and the total kinetic energy e\ they in- 
ject into the medium. These parameters and the 
OB star formation rate Sob completely charac- 
terize the forcing mechanism. It is important to 
emphasize that If is not the size ultimately reached 



by the expanding shells, but the scale over which 
the force is directly applied. The remainder of 
the shell expansion is inertial and strictly speak- 
ing docs not constitute a forcing. Our main con- 
clusions and implications are: 

1. The non-ubiquitous, small-scale nature of the 
forcing implies that the forced regions have small 
filling factors, giving rise to the coexistence of both 
forced and decaying turbulent regimes within the 
same flow. 

2. The kinetic energy injection and dissipation 
rates are always very close to each other, strongly 
suggesting that most of the dissipation occurs at 
or near the sources, where shocks are common, 
and that, for practical purposes, the injection and 
dissipation times are equal. In general, the flow 
dissipates its turbulent energy rapidly, in roughly 
15-20 Myr. 

3. In the forced regime, the global dissipation 
timescale t& is empirically found to depend mainly 
on the forcing scale If , and only very weakly on 
(or ek, the specific power injected per source) and 
Sob- 

4. Expressing t m in terms of average quantities in 
the flow and parameters of the forcing, we showed 
that id ~ u? ms /(ek /), where / is the filling fac- 
tor of the forcing regions. In terms of measurable 
properties of the ISM, t& > (S g )u 2 ms /(ek Sob), 
where (S g ) is the average gas surface density. For 
the solar neighborhood, > 1.5 x 10 7 yr. Since 
Sob oc S* oc S™, td oc Sg _1 (i?), td is expected to 
vary with Galactocentric radius. 

5. In the decaying regime, the turbulent kinetic 
energy decays as E^(t) oc (1 + t/i )~ a , with 
a w 0.8. For this value of a, the characteristic 
time for decay to half the initial energy is ~ 18 
Myr, but this time depends on the total kinetic 
energy content of the flow, because of the power- 
law dependence. We associate the decaying case 
to the regime existing in regions distant from the 
energy sources, which contains the "residual" tur- 
bulence left over from the strong local dissipation 
operating in the neighborhoods of the sources. 

6. From dimensional arguments, we suggest that 
the kinetic energy and the rms velocity disper- 
sion u lms of the "residual" turbulence (far from 
the forcing regions) should decay with distance 
£ as £-2«/(2-«) an d ^-«/(2-«) ) rcspcc tively. For 

a ~ 1, .Ek ~ i~ 2 and u lms ~ l~ x . 
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7. Our results imply that energy is dissipated 
much more efficiently at larger kinetic energy con- 
tents in the flow, in such a way that the dissipation 
time is rather insensitive to the total energy con- 
tent. 

8. If our results are applicable to the stratified 
vertical direction in the Galactic disk, then tur- 
bulent motions produced near the disk plane will 
propagate up to distances not too much smaller 
than the observed HI disk semi-thickness. This is 
consistent with models of galaxy evolution where 
the HI disk thickness is determined mainly by the 
turbulent kinetic energy content of the medium, 
which in turn results from the balance between ki- 
netic energy injection by stars and its dissipation 
rate. However, some extra support is probably 
necessary to keep the HI and ionized layers up at 
the observed heights. 

9. Our result of mostly local dissipation, again if 
applicable in the vertical direction, is in conflict 
with models of galaxy formation where the tur- 
bulent kinetic energy injected by SNe and stellar 
winds is assumed to reheat and drive back the gas 
from the disk into the cosmological dark matter 
halo (of sizes 15-20 times the disk size) in such a 
way that SF is self-regulated at the level of the 
whole intra-halo medium. Only in cases of non- 
stationary runaway SF (starbursts), most of the 
superbubbles might be able to blowout of the disk 
expelling large amounts of gas and energy into the 
intra-halo medium. 

7.2. Importance of small-scale, intermit- 
tent forcing 

The results of this paper imply that ISM tur- 
bulence is essentially different from ideal homo- 
geneous turbulence due to the nature of the in- 
tervening energy injection (forcing) mechanisms. 
Although several authors (e.g., Scalo 1987; Fleck 
1983; Norman & Ferrara 1996; see also the reviews 
in Vazquez-Semadeni 1999; Vazquez- Semadeni et 
al. 2000) have already discussed the fact that in 
"traditional" turbulence energy is injected at large 
scales while in the ISM it is injected over a range of 
scales (and preferentially at small ones), the con- 
sequences of the spatially- and temporally-discrete 
nature of interstellar energy injection have barely 
been discussed to our knowledge. 

The effects of intermittent forcing have recently 



begun to be considered in other frameworks as 
well. Bee, Frisch & Khanin (2000) have studied 
the statistics and solvability properties of Burg- 
ers turbulence in the presence of "kick" forcing 
(i.e., the forcing is a series of delta functions in 
time) , pointing out that such a flow combines fea- 
tures of both purely decaying and continuously 
forced cases. Leorat, Passot & Pouquet (1990) and 
Vazquez-Semadeni, Gazol & Scalo (2000) have 
shown that small-scale forcing is able to prevent 
the development of the gravitational and thermal 
instabilities, respectively. Finally, Kornreich & 
Scalo (2000) have considered the exposure of star- 
less interstellar clouds to the intermittent passage 
of shock waves originated at distant injection sites, 
to show that the typical time between passages 
should be comparable to the decay time within 
the clouds, explaining their apparently turbulent 
state even in the absence of local energy injection 
sources. 

In conclusion, the discrete, small-scale nature 
of interstellar kinetic energy injection seems to 
be responsible for a new, rich kind of turbulent 
flow, whose dissipation properties seem to pose 
clear constraints on models of galaxy formation 
and evolution. 

This work has benefitted extensively from com- 
ments from J. Brasseur, A. Brandenburg, E. Os- 
triker and J. Stone while one of us (E.V.-S.) at- 
tended the "Astrophysical Turbulence Program" 
of the Institute of Theoretical Physics of the Uni- 
versity of California at Santa Barbara, and from a 
thorough referee's report by M.-M. Mac Low. This 
work was supported by CONACYT grant 27752-E 
to E.V.-S. 
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Fig. 1. — Onc-dimcnsional profiles of the accelera- 
tion a x {x) {dotted line) and the final velocity after 
applying a x (x) over a stellar lifetime (AioB = 6.8 
Myr) for the parameters of the fiducial simulation 
labeled run 30. The amplitude a max is measured 
from peak to trough of the a x {x) curve. 
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Fig. 2. — Images of the density field (left) and the 
x-component of the acceleration for the fiducial 
run at time t = 272 Myr. The box size is 1 kpc and 
the resolution is 128 grid points per dimension. 



Fig. 3. — Evolution of various quantities for the 
fiducial run, run 30: total kinetic energy (dot- 
dashed line); total time derivative of the kinetic 
energy, dE^/dt (dotted line); kinetic energy injec- 
tion rate E, (solid line); kinetic energy dissipation 
rate E<j (dashed line), as defined by eq. (5), and 
the dissipation time id, defined by eq. (4) (dotted 
line with diamonds) . The latter is shown averaged 
over periods of 20 consecutive code outputs, due 
to its sensitivity on E i; and is given in units of 
10 8 years. The other quantities are given in code 
units. 
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Table 1 

Direct and indirect energy injection parameters 



Run # Resolution 81 Z f /pc Uf/km s" 1 e k b e k /10 49 crg c (S Q B) d / e 



6 (ISM) f 


128 


30.0 


- 7.0 






- 14.0 




12 (dec.) 8 


128 


NA 


NA 


NA 


NA 


NA 


NA 


17 (512 dec.) 


512 


NA 


NA 


NA 


NA 


NA 


NA 


30 (fiducial) 


128 


66.5 


7.02 


1.99 


2.84 


0.71 


.36 


33 


128 


38.0 


6.44 


2.82 


1.31 


0.49 


.10 


34 


128 


66.5 


3.21 


0.38 


0.54 


0.41 


.26 


35 


128 


66.5 


7.02 


1.99 


2.84 


1.63 


.55 


36 


128 


66.5 


3.21 


0.38 


0.54 


1.02 


.45 


37 


128 


38.0 


6.44 


2.82 


1.31 


1.13 


.18 


38 


128 


38.0 


3.21 


0.48 


0.23 


0.61 


.11 


39 


128 


38.0 


3.21 


0.48 


0.23 


0.28 


.06 


40 


128 


133. 


2.92 


0.22 


1.25 


0.58 


.67 


41 


128 


133. 


4.10 


0.69 


3.94 


0.71 


.68 


42 


128 


66.5 


7.02 


1.99 


2.84 


0.10 


.09 


43 


128 


38.0 


1.64 


0.12 


0.06 


0.45 


.09 


44 (512 fid.) 


512 


66.5 


5.86 


1.63 


2.33 


0.73 


.35 



a Number of grid points per dimension 

b In units of 10~ 3 erg s^ 1 gr _1 

c Total energy injected per stellar event 

d Time-averaged "OB-star" formation rate (in units of 10~ 5 yr _1 kpc~ 2 ), from t = 40.8 to t = 272 Myr. 
c Measured filling factor of all forcing regions 

f " — " denotes data not measured for this run 
S NA means not applicable 
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Table 2 

Global quantities measured in the runs 



rill 11 jf- 


/ 7Tiin\b 

\*\ / 




rms) /km S 1 


\£d// M yr 




iin(pred.) e 


A E /pc 




6 (ISM) h 


- 2.8 


1.80 


6.8 


18.0 


20.4 


NA 


245.' 


34.7 


12 (dec.) 


NA 


NA 


NA 


NA 


18.0 


NA 






17 (512 dec.) 


NA 


NA 


NA 


NA 


16.6 


NA 






30 (fiducial) 


2.36 


1.20 


5.5 


16.5 


16.3 


13.8 


163. 


27.4 


33 


0.80 


0.32 


3.0 


13.0 


12.8 


11.1 


107. 


34.4 


34 


0.33 


0.25 


2.6 


25.0 


24.3 


21.8 


178. 


66.1 


35 


5.43 


2.39 


7.6 


15.0 


14.2 


17.3 


151.5 


19.4 


36 


0.85 


0.50 


3.5 


19.0 


18.9 


23.3 


151.5 


45.0 


37 


1.56 


0.59 


4.0 


12.0 


12.1 


10.8 


107. 


26.2 


38 


0.15 


0.06 


1.3 


13.0 


13.0 


10.8 


108. 


84.5 


39 


0.07 


0.03 


0.9 


13.0 


13.2 


10.8 


107. 


115. 


40 


0.90 


0.86 


4.4 


33.0 


30.8 


48.4 


244.5 


55.1 


41 


1.84 


1.69 


6.2 


31.0 


29.6 


42.1 


271. 


46.4 


42 


0.40 


0.33 


3.0 


28.0 


26.4 


16.2 


173. 


50.3 


43 


0.03 


0.01 


0.6 


14.6 


14.0 


10.8 


104. 


178. 


44 (512 fid.) 


2.15 


1.07 


4.9 


16.0 


16.0 


13.7 


143.J 


26.6 



a All the measured quantities are averages over 40.8 < t < 272 (Myr) and correspond to a total area of 1 
kpc 2 

b Energy injection rate (in units of 10 36 erg/s) 
c Energy content (in units of erg) 

d Measured characteristic injection time scale (in Myr) 

e Predicted characteristic injection time scale according to eq. (10) (in Myr) 

f Length scale corresponding to the energy spectrum centroid (see eq. [20]) measured at t = 136 Myr (51.7 
Myr for run 44) 

s Large-scale crossing time t CI = A e /w rms (in Myr) 

h " — " and NA mean not measured and not applicable, respectively 

'Centroid calculated at t = 272 Myr to allow for a roughly stationary regime in this run 

J Centroid calculated at t — 68.0 Myr (maximum evolution time for this run) 
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Fig. 4. — Scaling of the rms velocity dispersion 
Wrms vs. Ek for run 6. The points define a line 
with a slope very close to 1/2. 




Fig. 5. — Images of the log of the density field of 
run 6, which used density-directed star formation 
(SF), self-gravity and the Coriolis force, shown at 
t = 15.0 Myr (left) and at t = 149.6 Myr (right). 
SF started at t = 5.2 Myr in this run. In the 
left panel small expanding shells are seen, still in 
their initial phases. The trapezoidal void near the 
center was not formed by stellar activity, but by 
the turbulent initial conditions. At the later time, 
large shells of up to ~ 500 pc are seen. These are 
due to induced SF in the shells. This run exhibits 
a similar dissipation time as the fiducial run (30) 
in spite of the more realistic conditions. 




Fig. 6. — Log of the density field for the decay sim- 
ulation, run 12, at times t = (left) and t — 189.0 
Myr (right). This runs is a restart of run similar 
to run 6 at t = 49.0 Myr, with SF, the Corio- 
lis force and self-gravity turned off. Expanding 
shells can still be seen at the earlier time, while 
a much smoother density structure is seen at the 
later time. 
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0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 
log (1+t/t.) 

Fig. 7. — Evolution of the total kinetic energy (in 
code units) for the decaying run (run 12). 
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Fig. 8. — Evolution, for run 17 (similar to run 12 
but at resolution 512 2 ), of the log of the kinetic en- 
ergy (solid line with plus signs, in code units) and 
of the absolute value of the most negative diver- 
gence of the velocity field (dashed line, displaced 
by —2 in the vertical direction to fit it in the plot). 
The dash-dotted line shows a slope of —0.9. 




Fig. 9. — a) (Left) Velocity spectra of the fiducial 
run at times t — 1.36 Myr lower line, t — 6.8 Myr 
(middle line) and t = 136 Myr (upper line), b) 
(Right) rr-component of the velocity field at t = 6.8 
Myr. At t = 1.36 Myr, the spectrum is domi- 
nated by the (mild) initial conditions, but a peak 
at logfc = 1.4 corresponds to the first (and only) 
shell present at that time, just formed. At t = 6.8 
Myr this peak has moved to log k ~ 1.2 as the shell 
has expanded, and several other ripples at larger k 
correspond to new shells. At lower k the spectrum 
corresponds to the ensemble of shells. At t = 136 
Myr the spectrum has reached a power-law shape, 
indicating fully developed turbulence. 
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